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ABSTRACT 

Extremely large opaque troughs in the Lya forest have been interpreted as a sign of an ex¬ 
tended reionization process below z ^ 6. Such features are impossible to reproduce with 
simple models of the intergalactic ionizing background that assume a uniform mean free path 
of ionizing photons. We build a self-consistent model of the ionizing background that includes 
fluctuations in the mean free path due to the varying strength of the ionizing background and 
large-scale density field. The dominant effect is the suppression of the ionizing background in 
large-scale voids due to “self-shielding” by an enhanced number of optically thick absorbers. 
Our model results in a distribution of 50 Mpc//i Lya forest effective optical depths that sig¬ 
nificantly improves agreement with the observations at z ~ 5.6. Extrapolation to z ~ 5.4 and 
z ~ 5.8 appears promising, but matching the mean background evolution requires evolution in 
the absorber population beyond the scope of the present model. We also demonstrate the need 
for extremely large volumes (> 400 Mpc on a side) to accurately determine the incidence of 
rare large-scale features in the Lya forest. 


1 INTRODUCTION 

The reionization of hydrogen was an important milestone in the his¬ 
tory of the Universe, representing the culmination of early structure 
formation and a dramatic phase change in the intergalactic medium 
(IGM). Great efforts have been made to investigate the reionization 
epoch, both theoretically and observationally, as it provides a pow¬ 
erful tool for testing theories of the formation of the first galaxies 
(Loeb & Furlanetto 2013). 

One of the dominant features of neutral hydrogen is strong 
Lyman-series absorption. Observations of large-scale opaque re¬ 
gions in the Lya forest of high-redshift quasars, also known as 
Gunn-Peterson troughs (Gunn & Peterson 1965), have already 
placed intriguing constraints on the end stages of reionization (e.g. 
Fan et al. 2001, 2006; Mesinger 2010; McGreer et al. 2011, 2015). 
Searches for Lya damping wing absorption due to the presence of a 
neutral IGM (Miralda-Escude 1998) in high-redshift quasar spectra 
have been inconclusive (e.g. Mesinger & Haiman 2007; Mesinger 
& Furlanetto 2008b; Schroeder et al. 2013; Bolton et al. 2011; Sim- 
coe et al. 2012 but see Bosman & Becker 2015), but could eventu¬ 
ally provide “smoking gun” evidence for primordial neutral ma¬ 
terial. Lya damping wing absorption due to a neutral IGM may 
also be responsible for the precipitous drop in the fraction of star¬ 
forming galaxies that show strong Lya emission lines above z ~ 6 
(e.g. Stark et al. 2010; Ono et al. 2012; Schenker et al. 2014; Pen- 
tericci et al. 2014; Tilvi et al. 2014), subject to considerable model- 
dependent uncertainties due to the complexity of Lya emission and 
scattering processes and inhomogeneous reionization (e.g. Santos 
2004; Furlanetto et al. 2004; McQuinn et al. 2007; Mesinger & 
Furlanetto 2008a; Dijkstra et al. 2011; Bolton & Haehnelt 2013; 
Choudhury et al. 2014; Mesinger et al. 2015). 

The ionization state of gas in the IGM is regulated by the “ion¬ 


izing background”, the extragalactic ionizing radiation field. The 
ionizing background is comprised of ionizing photons emitted by 
young stars and quasars, filtered by neutral gas structures that set 
an average attenuation length or mean free path A (e.g. Haardt & 
Madau 1996, 2012; Faucher-Giguere et al. 2009). Interpretation of 
Lya forest observations typically assumes a uniform hydrogen ion¬ 
ization rate Fh i due to the large mean free path (A > 100 Mpc) 
of ionizing photons at z < 5 that smooths the radiation field on 
large scales. More sophisticated models include the effect of dis¬ 
crete sources of ionizing photons and attenuate ionizing radiation 
with a constant mean free path (e.g. Meiksin & White 2004; Mc¬ 
Donald et al. 2005; Wyithe & Loeb 2006; Mesinger & Furlanetto 
2009). Recently, Becker et al. (2015) (henceforth B15) discovered 
an enormous ~ 100 Mpc/fi opaque trough in the Lya forest spec¬ 
trum of ULAS J0148-I-0600 covering z ~ 5.5 —5.9 which is at odds 
with their uniform mean free path ionizing background simulations. 
They suggest that this discrepancy could be due to fluctuations in 
the mean free path of ionizing photons, possibly as a result of an 
extended reionization process. 

However, the mean free path of ionizing photons is likely to 
vary substantially at z ~ 5.6 even if ionization equilibrium is as¬ 
sumed. Post-processing of hydrodynamical simulations of the IGM 
suggests that the mean free path should vary with the strength of 
the ionizing background as A oc F^j with ^ ~ 0.6-0.8 (Mc¬ 
Quinn et al. 2011; henceforth Mil), and uniform mean free path 
ionizing background simulations at z ~ 5.6 lead to large-scale 
(tens of Mpc) regions with enhanced and diminished background 
(Mesinger & Furlanetto 2009). Regions with a weak ionizing back¬ 
ground should thus have a short mean free path, while regions with 
a strong ionizing background should be more transparent to ioniz¬ 
ing photons. This effect causes a feedback loop that should lead to 
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enhanced contrast of the ionizing background on (roughly speak¬ 
ing) scales larger than the average mean free path of ionizing pho¬ 
tons. Fluctuations in the mean free path have recently been studied 
with linear theory at 2 ~ 2.5 (Pontzen 2014; Gontcho A Gontcho 
et al. 2014) in terms of subtle modulation of the baryon acoustic 
oscillations signal, but have been ignored at higher redshift where 
they are potentially very important to the large-scale structure of 
the radiation field. 

In this work, we construct a semi-numerical model to compute 
the ionizing background in a large cosmological volume at 2 ~ 5.6 
including a varying mean free path of ionizing photons. The local 
mean free path is iteratively determined by the local overdensity 
and strength of the ionizing background, the latter of which is com¬ 
puted by filtering ionizing radiation through the varying opacity 
field. The resulting radiation field is “self-consistent” in the sense 
that the opacity of the IGM and ionizing background are set by 
each other - the opacity depends on the strength of the ionizing 
background, and the ionizing background itself is filtered by that 
opacity. We then investigate the statistics of large-scale features in 
Lya forest transmission that our new ionizing background model 
implies and compare to the observations compiled by B15. 

The structure of the paper is as follows. In Section 2, we de¬ 
scribe our models for sources and absorbers of ionizing photons. In 
Section 3, we compute the self-consistent fluctuations in the ion¬ 
izing background and mean free path. In Section 4, we apply our 
ionizing background model to the statistics of large-scale Lya for¬ 
est absorption. In Section 5, we test the redshift evolution of our 
model across 2 ~ 5.8-5.4. In Section 6, we discuss the implica¬ 
tions of our model for the IGM and sources of ionizing photons at 
high-redshift. Finally, in Section 7 we conclude with a summary 
and speculate on future observational tests of our model. 

In this work we assume a standard ACDM model with Gm = 
0.3, Ga = 0.7, h = 0.7, and erg = 0.82, consistent with con¬ 
straints from the cosmic microwave background (Hinshaw et al. 
2013, and close to the latest bounds by Planck Collaboration et al. 
2015). Distance units are comoving unless otherwise noted. 


2 SOURCES AND SINKS OF IONIZING PHOTONS 

Our three-dimensional model for the ionizing background requires 
two ingredients: the distribution of ionizing sources and a prescrip¬ 
tion for spatially variable IGM opacity to ionizing photons (i.e. a 
varying mean free path). 

2.1 Semi-numerical Density and Halo Fields with DEXM 

To construct the distribution of ionizing sources in our model, 
we use the semi-numerical cosmological simulation code DEXM 
(Mesinger & Furlanetto 2007; henceforth MF07). In the following, 
we summarize the method employed in DEXM, deferring a dis¬ 
cussion of the detailed implementation to MF07. We begin with 
a 2100^ grid of Gaussian cosmological initial conditions (ICs) in 
a volume 400 Mpe on a side. While this volume is considerably 
larger than the majority of cosmological simulations, we will in¬ 
vestigate whether this is a representative volume in Section 6.2. 
The linear velocity field is then computed on a coarser grid of 
700^ cells using the ZeFdovich approximation (ZeTdovich 1970; 
Efstathiou et al. 1985). In the public release of DEXM, the quasi- 
linear density field is then computed by displacing “particles,” cor¬ 
responding to the grid of initial conditions, following the coarsely- 
computed velocity field (with no interpolation) and assigning their 


mass to the nearest grid cell to compute the density field. The right 
panel of Figure 1 shows that while the standard DEXM density field 
is reasonably smooth in high-density regions with a large number 
of particles per cell, low-density regions suffer from considerable 
shot noise. In particle-based cosmological simulations, a nearest- 
neighbor smoothing approach is typically performed on the parti¬ 
cle distribution to estimate the continuous density field (e.g. Mon¬ 
aghan 1992; Springel 2005). Nearest-neighbor schemes naturally 
degrade resolution in low-density environments in favor of high- 
density regions, ideal for computationally intensive hydrodynami- 
cal simulations of galaxy formation. However, these void environ¬ 
ments are crucial to modeling transmission through the Lya forest 
at 2 > 5 (Bolton & Becker 2009), so we adopt a novel approach 
that does not degrade our simulation resolution through explicit 
spatial smoothing. 

The linear velocity field produced by the ZeFdovich approxi¬ 
mation is typically smooth on small scales, its structure dominated 
by large-scale bulk flows of material into sheets and filaments. To 
compute the quasi-linear density field, we displace a super-resolved 
grid of 21000^ IC particles^ by interpolating the ZeFdovich ve¬ 
locity field via trilinear interpolation and then binning these sub¬ 
particles onto a 700^ grid.^ The material in each original IC grid 
cell is not only displaced but also stretched and sheared due to gra¬ 
dients in the velocity field. Because the velocity field is smooth 
on small scales, this procedure results in a smooth distribution of 
matter on a uniform grid without any additional spatial filtering - 
even in the lowest density environments. The left panel of Figure 1 
shows the resulting density field computed from the same ICs as 
the right panel, demonstrating the lack of artifacts - especially in 
void environments - compared to the standard method. 

We use the halo finding algorithm of DEXM to populate the 
simulation volume with a realistically-clustered distribution of dark 
matter halos. In brief, halos are located by filtering the linear den¬ 
sity field on successively smaller scales and searching for (non¬ 
overlapping) regions where the filtered density 5{x, M) is greater 
than the linear collapse threshold 5c- We locate halos with masses 
as low as Mf, ~ 2 x 10® Mq (~ 30 IC cells) and show the result¬ 
ing mass function at 2 = 5.6, which (by construction) agrees very 
well with the Sheth & Tormen (1999) mass function in Figure 2. We 
then displace the halos using the ZeFdovich velocity field described 
above, enhancing their bias on large-scales and resulting in cluster¬ 
ing statistics that are consistent with N-body simulations (MF07). 
As discussed by MF07, we do not expect one-to-one agreement be¬ 
tween halo locations found with this procedure and those formed 
in an N-body simulation with the same initial conditions, but the 
agreement with large-scale clustering statistics is adequate for our 
purposes. 

To compute the ionizing photon output of each halo, we abun¬ 
dance match (e.g. Vale & Ostriker 2004) to the observed UV lu¬ 
minosity functions from Bouwens et al. (2015b), interpolating be¬ 
tween their 2 ~ 5 and 2 ~ 6 best-fit Schecter functions and as¬ 
suming a constant ratio of /ion/j49i2 between the non-ionizing and 
ionizing continuum luminosity as in B15, where A 912 = 6.0 is 

^ The required number of sub-pailicles to produce a qualitatively smooth 
(and converged) final density field depends on the relative number of IC vs. 
velocity grid cells, and on the spatial resolution of the output density field. 

^ Matching the resolution of the ZeFdovich velocity field is not strictly 
required - one could in principle bin the sub-particles onto a coarser or 
finer grid. In the latter case, the resulting density field appears to have subtle 
box-like artifacts, presumably due to the linear interpolation of the velocity 
field. 
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Figure 1. Comparison of the standard nearest-gridpoint Zel’dovich approximation density field from DEXM (left) with our interpolated Zel’dovich velocity 
approach (right). The slice shown is 0.57 Mpc thick (1 pixel) and 100 Mpc on a side. 



Figure 2. Sheth & Tormen (1999) halo mass function at 2 = 5.6 (solid 
curve) compared to the mass function of halos in our DEXM simulation (open 
points). 


the expected ratio for a young stellar population and /ion is a free 
parameter that represents a combination of the escape fraction of 
ionizing photons /eac and uncertainties in the value of ^ 4912 . For 
comparison to the B15 models we further assume that the spec¬ 
trum of ionizing photons emitted by each galaxy is a power law 
Ln oc with a = 2.0.^ 


® See Section 5.1 of Becker & Bolton (2013) for a detailed assessment of 
this choice. 


2.2 Varying Mean Free Path 

The mean free path of ionizing photons in the IGM has been mea¬ 
sured out to 2 ~ 5 through stacking analyses of quasar ioniz¬ 
ing continua (Worseck et al. 2014) and by counting absorption 
lines in the Lya forest (e.g. Songaila & Cowie 2010; Rudie et al. 
2013). Typically this mean free path is assumed to be uniform in 
space. Flowever, the mean free path itself should be sensitive to 
the strength of the ionizing background due to the regulation of 
Ft I absorber sizes (Mil; see also Munoz et al. 2014) and overdense 
regions likely contain more dense self-shielded gas. 

Mil developed an analytic framework to describe the self- 
consistent relationship between the mean free path A and the ion¬ 
ization rate F based on the IGM model of Miralda-Escude et al. 
(2000) (henceforth MHR). The MHR model describes the ioniza¬ 
tion state of the universe by assuming that all gas with density 
above a critical value Ai is neutral, self-shielded from the ioniz¬ 
ing background, while all lower density material is completely ion¬ 
ized.^ This means that global parameters that depend on the amount 
of neutral or ionized gas can be completely specified through the 
gas density probability distribution function (PDF) f’(A). Ad¬ 
ditionally, the MHR model assumes that the mean free path is 
given by A oc Fy(Ai)”^'^®, where Fv{Ai) is the volume fill¬ 
ing factor of gas above the density Ay If the gas density PDF 
can be approximated by a power law P{A) oc A“^ in the den¬ 
sity regime that dominates the IGM opacity, then A oc 
Under the assumption of ionization equilibrium. Mil showed that 
F oc A^^ In other words, A oc Fundamen¬ 

tally, if the ionization rate is larger, higher density gas will be ion¬ 
ized, reducing the total volume of neutral gas and thus increasing 
the mean free path (see also Munoz et al. 2014). Mil then con¬ 
firmed this scaling by performing radiative transfer through hydro- 
dynamic cosmological simulations. 

In detail, a substantial amount of IGM ionizing opacity should come from 
optically thin absorbers. However, the scaling of A with F should not be 
sensitive to this assumption, as confirmed in numerical simulations by Ml 1. 
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In our fiducial simulations we allow the mean free path to vary 
as A oc corresponding to a gas density PDF -P(A) oc 
This PDF is equivalent to assuming that all gas lies in spherically 
symmetric clouds with an isothermal density profile, i.e. p oc r~^. 
This is likely a conservative assumption; Mil found that the true 
PDF is likely to be steeper at high redshift, leading to a more sen¬ 
sitive dependence with Fhi. In Section 6 we will discuss the effect 
of modifying this choice. 

The variation of the mean free path with large-scale density 
depends on the relative number density of self-shielded absorbers 
as a function of that large-scale density. In general, this variation 
should moderate the amplitude of mean free path fluctuations be¬ 
cause (if absorbers are biased relative to dark matter) it acts in the 
opposite direction from the F dependence described above (in that 
voids will become more transparent). At z > 5.5, the overdensity 
threshold for self-shielding is less than the virial overdensity, and 
so the population of absorbing clouds is likely dominated by the 
distant outskirts of very low-mass halos and filamentary structures 
in the IGM (Munoz et al. 2014). These overdense structures are un¬ 
resolved in our simulations and so we only attempt to describe the 
variation in their abundance on the same coarse scale as the ioniz¬ 
ing background (see below). The absorbing clouds are likely only 
weakly biased with respect to the density field, so we assume a fidu¬ 
cial absorption bias of unity for simplicity (i.e. A oc A“^). How¬ 
ever, in the densest environments in our simulations which have 
the strongest ionizing background, the absorber population may in¬ 
stead be dominated by gas within the virial radius of low-mass ha¬ 
los Munoz et al. (2014), leading us to overestimate the mean free 
path in these regions. We leave a full exploration of the relationship 
between large-scale overdensity and the mean free path for future 
work. 

For comparison with B15, we assume that the column den¬ 
sity distribution of absorbers that are most relevant to H l opacity 
(AhI ~ All ~ cm~^) behaves as a power law 

/(Ah i) oc A“'® with /3 = 1.3, consistent with the 2 ~ 2-5 model 
in Becker & Bolton (2013). This implies a frequency dependence of 
the mean free path of A^ oc oc ® (e.g. Haardt & Madau 

1996). In practice, our results are insensitive to this choice because 
of the already steep frequency dependence of the H l cross section 
(oc u~^) and galaxy spectrum (oc i'~^) which also contribute to the 
ionization rate calculation in the following section. 

The final expression that we use for the mean free path is then: 


A(r, A, p) = Ao(rHi/ro)"/"A-i(i//t.Hi)°®, (i) 


where Ao and Fo are tuned, via re-normalizing the ionizing emis- 
sivity with /ion, such that an ionizing background model with a 
spatially constant A = Aq produces (Fhi) = Fq. 

Extrapolation of the Worseck et al. (2014b) mean free path 
evolution fit (derived from measurements over 2 ~ 2.5-5.2) to 
2 ~ 5.6 suggests A ~ 54 Mpc as an estimate of the average mean 
free path. However, the ionization rate in the IGM is constant from 
2 < 2 < 5 (Becker & Bolton 2013) before dropping by roughly an 
order of magnitude from 2 ~ 5 to 2 ~ 6 (Bolton & Haehnelt 2007; 
Wyithe & Bolton 2011; Calverley et al. 2011). This suggests that 
extrapolation of this fit to higher redshifts is likely to overestimate 
the mean free path by as much as a factor of a few. In this work 
we are interested in exploring a range of fluctuating background 
scenarios, so we investigate models with Ao = 15, 22, and 34 Mpc. 


3 NUMERICAL MODEL OF THE IONIZING 
BACKGROUND 


In this section, we describe the method used to compute the ioniz¬ 
ing background given the ingredients listed in the previous section. 

Our fiducial ionizing background models are computed on 
a coarse 80® grid, corresponding to cells of 5 Mpc on a side. 
This coarse resolution is sufficient for our purposes because it re¬ 
solves the typical mean free path in our simulations (which sets the 
smoothing scale) and we have found little difference in the large- 
scale radiation field structures over a range of grid sizes from 48® 
to 100®. We also bin the density field onto the same coarse grid 
for computational efficiency - as described in the previous sec¬ 
tion, we do not resolve neutral gas in our simulations, but instead 
assume that the number density of absorbers scales proportional 
to the large-scale overdensity. The halo-finding procedure in Sec¬ 
tion 2.1 results in ~ 2.4 x 10^ halos. While a uniform mean free 
path calculation of the radiation field is computationally trivial for 
any number of sources (i.e. J,^ oc ^ /R^, or computed 

via FFT), self-consistently including the mean free path variations 
from Section 2.2 requires line-of-sight integrations between each 
source and grid cell. For computational efficiency, we bin the ion¬ 
izing sources from Section 2.1 - originally found on a 2100® grid 
- onto the ionizing background grid, then compute the local and 
non-local contributions to the ionizing background separately.® 

The non-local contribution to the specific intensity J„ at cell i 
is summed over every other cell j located at fj , 


Ju.i — 




^ (47r|ri -rj|)2 


( 2 ) 


where is the specific luminosity of cell j and T„{ri, fj) is the 
integrated optical depth between cells i and j. 


T^{ri,fj)= / [Ai,(rHi(r), A(r))]"^dr, 


(3) 


where Aj/ is the mean free path of frequency v photons correspond¬ 
ing to Fhi and A along the line of sight. By integrating along the 
line of sight, the varying H l opacity of the IGM due to fluctuations 
in the strength of the radiation field and the large-scale density field 
is explicitly taken into account as described in the previous sec¬ 
tion.® 

The local contribution is computed by assuming that sources 
are spread evenly throughout each cell and can be described with a 
uniform emissivity e^. The radiation field at the center of a uniform 
emissivity sphere of radius I /2 and mean free path A is simply the 
integral over spherical shells with luminosity dL = e x Trrr^dr 
or 


J local,sphere 
u 



e^(47rr®) 

(47rr)2 


dr 


47r 


(4) 

which reverts to the “absorption limited” approximation for the ion¬ 
izing background intensity J eA/dvr (Meiksin & White 2003) 


^ This reduces computation time in two ways: the number of sources drops 
from ~ 2.4 x 10^ to ~ 5 x 10^, and the opacity between each pair of cells 
only has to be computed once. Thus, the total computation time is decreased 
by a factor of ~ 100. The gain in efficiency depends on the resolution of 
the 3D grid. 

^ In this simple model we neglect cosmological effects, most notably the 
redshifting of ionizing photons as they travel from the source cell. Fortu¬ 
nately, the assumed mean free path of ionizing photons is short enough that 
this is unlikely to have a significant effect on our results. 
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when integrated out to infinity. There is no simple exact expression 
for the radiation field at the center of a cubical volume with uniform 
emissivity, but we find numerically that the approximation 


J local,cube 
u 


Ttt 


(5) 


with C, ~ 0.72 results in an average ionization rate that is indepen¬ 
dent of spatial resolution. 

Finally, the hydrogen ionization rate Fh i in each cell is com¬ 
puted by integrating over frequency, 


FHI.i 



7 ■ -I- 7‘‘ 
hv 




( 6 ) 


where cthi(i') is the hydrogen photoionization cross-section from 
Verner et al. (1996). This additional integration step considerably 
increases the computation time of our simulations relative to uni¬ 
form mean free path models. By testing with lower-resolution 48^ 
models, we find that the ionizing background fluctuations in the 
frequency-integrated model can be well-reproduced by a single¬ 
frequency model with v « 1.32tzHi, and use this approximation 
to speed up computation of our highest resolution models.^ 

Because the local mean free path is a function of Thi, the 
above procedure must be iterated ~ 10 times until the average 
change of Th i per cell is less than 0.3%. At this point the ioniza¬ 
tion rate in the vast majority of the volume has converged, except 
for a very small fraction of cells at the center of opaque void regions 
that are strongly shielded from all ionizing sources. 

In the rest of the paper, we will describe a series of mod¬ 
els with a varying average mean free path. Models that employ a 
uniform mean free path (i.e. = constant in equations 2 and 3 

above) will be referred to as “uniform-A” while those that include 
our full prescription for mean free path variations will be referred to 
as “fluctuating-A”. The uniform-A models have been normalized to 
reproduce a single average F by tuning the normalization of the ion¬ 
izing emissivity through /ion (analogous to changing the assumed 
escape fraction of ionizing photons), and the fluctuating-A mod¬ 
els were computed using the uniform-A model as the first iteration. 
Throughout most of the paper, we will utilize a single realization of 
the density field (with a volume of [400 Mpc]®) in order to illumi¬ 
nate the dependence of the results on the absorber properties. We 
will examine the role of cosmic variance in Section 6.2 separately, 
but we note that our fiducial realization appears to be a conservative 
estimate of the impact of a fluctuating mean free path. 


3.1 Ionizing Background Results 

In Figure 3, we show the distribution of Fhi for our fiducial 
fluctuating-A models assuming A = 15, 22, and 34 Mpc. The solid 
curves show the self-consistent fluctuating-A models, while the dot¬ 
ted curves show the corresponding uniform-A models. Mean free 
path fluctuations substantially broaden the distribution of Fh i to¬ 
wards lower values, and the strength of the effect strongly depends 
on the average A. This tail to low Fhi comes about in large-scale 
underdensities where the weak Fhi - and thus, the shorter A - acts 
to “self-shield” the region from distant sources of ionizing photons 
which would otherwise lead to a more uniform, shallower trough 


^ The 80® single-frequency models presented here require approximately 7 
hours of runtime per iteration on a 12-core Intel Xeon (Nehalem generation) 
computer. 



Figure 3. The solid curves show FHi/jFHi) from the fiducial ionizing 
background simulation at 2 : = 5.6 for A = 15, 22, 34 Mpc in red, blue, 
and orange, respectively. The dotted curves show the distributions from the 
corresponding uniform MFP models (i.e. the first iteration of the ionizing 
background calculation). The differences between the two sets therefore 
demonstrate the effect of mean free path fluctuations. 


in Fhi- The high-Fni end of the distribution is only modestly ad¬ 
justed for two reasons. First, the density and ionization adjustments 
to the mean free path partially cancel out, leading to only a modest 
boost in the mean free path. Second, close to galaxies the structure 
of the radiation field is dominated by the 1 /r® term in Ji^ (equa¬ 
tion 2), so increasing the mean free path has a lesser effect than one 
might expect (i.e. oc A, see Meiksin & White 2003). 

Figure 4 shows density, halo, and ionizing background fields 
centered on the one slice of our fiducial model. The bottom panels 
demonstrate the spatial coherence of the background fluctuations 
with 5 Mpc (1 pixel) thick slices of the A = 15 Mpc model in 
the uniform-A (left) and fluctuating-A (right) cases. The addition 
of a fluctuating mean free path greatly increases the contrast of the 
radiation field on > A scales. The upper panels of Figure 4 show 
projected density and halo fields ±10 Mpc from the center of the 
ionizing background slice. Comparison to the bottom panels shows 
that the spatial fluctuations in our ionizing background model are 
due to the inhomogeneous distribution of ionizing sources related 
to the formation of large-scale structure. Figure 5 shows the fluc¬ 
tuating mean free path in a slice corresponding to the lower right 
panel of Figure 4. In this case, A represents the “effective” mean 
free path inside each UVB pixel, i.e. the optical depth to hydrogen 
ionizing photons across the pixel is thi = dR/X where dR is the 
path length through the pixel. Large-scale features in the radiation 
field dominate the structure of mean free path map, with small-scale 
“noise” reflecting variations in the density field on the 5 Mpc pixel 
scale. 
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Figure 4. Top left: 20 Mpc-thick slice of the density field in our fiducial cosmological volume, 400 Mpc on a side. Top right: Halos found using the dexm 
halo-finding procedure inside the same slice. The size and color of each point represents the corresponding UV magnitude of the halo when matched to the 
Bouwens et al. (2015b) luminosity function, from Muv ~ —18 (black, smallest) through Muv ~ —22 (red, largest) in steps of AM = 1. The halos shown 
represent the brightest ~ 1% of all galaxies in this slice, with halo masses ranging from ~ 10^^'^-10 ^^'^Mq. Bottom left: Fluctuations in the uniform- 
A = 15 Mpc ionizing background model in a 5 Mpc-thick slice corresponding to the center of the slice shown in the density and halo fields. Variations of a 
factor of ~ 2-3 are common on large scales, as found by previous authors. Bottom right: Fluctuations in the fluctuating-A = 15 Mpc ionizing background 
model in the same 5 Mpc-thick slice. The addition of mean free path fluctuations greatly increases the fluctuations in Fhi, especially to very low values in 
large-scale underdensities. 


4 IMPLICATIONS FOR LARGE-SCALE Lya FOREST 
TRANSMISSION 


Our simulations do not resolve the density and velocity fields at 
the Jeans scale of the gas at these redshifts, so we cannot pro¬ 
duce realistic Lya forest spectra. We can still estimate the effect the 
ionizing background fluctuations on the large-scale opacity of the 
Lya forest by employing the fluctuating Gunn-Peterson approxi¬ 
mation (FGPA; e.g. Weinberg et al. 1997) on sightlines through our 
evolved density field, 


TOP 


■ VTSOOkJ is X 10-13 s-ij 

A 2-0.724(7-1) / I + 

V 6.6 ) 


(7) 


where To is the temperature of the IGM at the mean density, 7 is 
the slope of the temperature-density relation, and k is a normaliza¬ 


tion factor (e.g. Dixon & Furlanetto 2009) that we adjust to match 
the observed distribution of Lya forest optical depths, similar to 
B15®. We assume To = 7500 K and 7 = 1.5, consistent with 
standard models for the thermal history of the IGM® (e.g. Puch- 
wein et al. 2015). We choose a characteristic scale of 50 Mpc//i 
for comparison to B15 and compute the effective optical depth 
Teff = —ln(^ in steps of ~ 0.06 Mpc along each 

sightline. Our fiducial simulations have {Fh!) ~ 3 x 10“^® s”! 


3 In detail, B15 instead varied Fh i in their simulations to match observa¬ 
tions, but in practice the behavior is identical. 

® Such models typically assume that reionization is complete by 2 : > 9. If 
instead reionization did not end until 2 ~ 6, the thermal state of the IGM 
could be very different. The principal effect would be to increase the mean 
temperature and hence shift the required normalization factor k closer to 
unity. 
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Figure 5. Mean free path of hydrogen-ionizing photons in the ionizing 
background model slice shown in the bottom right panel of Figure 4, where 
A oc and A((rHi)) = 15 Mpc. The strong fluctuations in Fhi 

dominate the large-scale features, while small-scale “noise” is due to the 
relatively small variations in A. 

and require k ~ 0.14 to match the P(refr) observations by B15. 
The need for such a small value of k is primarily due to the inabil¬ 
ity of our simulations to resolve the low-density gas in small-scale 
cosmic voids which dominates the transmission in the relatively 
opaque z > 5 Lya forest (Bolton & Becker 2009). For the rest of 
this work, we assume that the dependence of Teff on Fhi is cap¬ 
tured by our relatively low resolution simulation and note that we 
cannot place any constraints on the true mean value of Fh i. 

In Figure 6 we show “maps” of Lya forest opacity along 50 
Mpc /h sightlines perpendicular to the page, centered on the same 
slice shown in Figure 4, for three different models: uniform F, 
uniform-A=15 Mpc, and fluctuating-A=15 Mpc from left to right. 
It is clear that coherent large-scale structure in the radiation field 
at this redshift plays an important role in the large-scale opacity of 
the Lya forest, and that the addition of mean free path fluctuations 
greatly enhances this effect - especially on scales larger than the 
average mean free path itself. 

The uniform Fhi slice in the left panel of Figure 6 is sim¬ 
ply a reflection of the projected density field - regions with higher 
density show less transmission, and vice versa. This is the standard 
picture of the Lya forest at lower redshifts when the ionizing back¬ 
ground is largely uniform. High density regions only make up a 
small portion of the volume in the standard cosmological model, so 
opaque sightlines are rare. Once the correlation between the density 
field and the radiation field is strong enough, as in the fluctuating-A 
model in the right panel of Figure 6, this picture reverses: low den¬ 
sity regions become opaque due to a dearth of ionizing photons, at 
least on scales larger than the average mean free path. 

To determine the statistical properties of large-scale features 
in the Lya forest throughout our simulation volume, we computed 
the Lya effective optical depth along 250000 randomly positioned 
and oriented 50 Mpc /h sightlines in our fiducial realization. In Fig¬ 
ure 7, we compare the observations by BI5 to our A = 15 Mpc 
models with and without mean free path fluctuations. The uniform 
background and uniform-A models are nearly identical, despite the 
clear difference in physical environments corresponding to a given 
optical depth seen in Figure 6. In agreement with B15 we find that a 
uniform ionizing background is sufficient to describe the observed 


P{< Tefs) < 0.5, but it fails to explain the tail to higher optical 
depths. Fluctuations in the mean free path naturally lead to a tail of 
high optical depths that are more consistent with observations than 
models that assume a uniform mean free path. 

In Figure 8 we show the distribution of Teir for A = 15,22, and 
34 Mpc fluctuating-A models. The highest Teir observed by B15 is 
extremely rare unless the average A is small, and it is still quite rare 
even in our A = 15 Mpc model (but see Section 6.2). We also note 
that the A = 15 Mpc model shows some subtle deviations from ob¬ 
servations at small TeS- These deviations come from regions expe¬ 
riencing a strong radiation field, i.e. environments close to galaxies, 
which our simplified density field model is unlikely to model accu¬ 
rately. For this reason, our model cannot be rigorously and quanti¬ 
tatively compared to the observations. Such an effort would require 
an improved treatment of the small-scale environments of the ioniz¬ 
ing sources, in addition to our large-volume treatment that matches 
the observed tail to high Teir. 

In Figure 9 we compare the relationship between the average 
density along each sightline and the resulting optical depth for the 
uniform ionizing background, uniform-A, and fluctuating-A mod¬ 
els. While the average density and optical depth are tightly corre¬ 
lated when the ionizing background is uniform, the situation be¬ 
comes more complex when the radiation field strongly fluctuates. 
Although the distribution of Tefi is very similar between the uni¬ 
form background and uniform-A cases (Figure 7), the physical na¬ 
ture of a sightline with a given Teir is different: in the uniform-A 
model the rare sightlines with high optical depth are modestly un- 
derdense and the densest regions are typically close to the median 
optical depth. This difference becomes greatly enhanced once mean 
free path fluctuations are included. 

We note that the distribution of TeS in our uniform background 
model (i.e. the left panel of Figure 6 and the solid blue curve in Fig¬ 
ure 8) is somewhat broader than in previous works (Becker et al. 
2015; Chardin et al. 2015), likely due to our poor physical reso¬ 
lution and relatively large cosmological volume (see the Appendix 
of Becker et al. 2015). Note, however, that the broader distribution 
seen in our simulations makes our estimate of the impact of a fluc¬ 
tuating mean free path a conservative one. Recall that the introduc¬ 
tion of a mean free path inverts the relation between density and 
optical depth, which is a positive correlation in the spatially con¬ 
stant F case (i.e., the blue contours in Figure 9). If that relation is 
weaker for a spatially constant ionizing background, the introduc¬ 
tion of the mean free path will result in even stronger fluctuations. 
Testing this in detail is beyond the scope of this paper, as it requires 
numerical simulations with extremely high dynamic range. 


5 EVOLUTION OF THE EFFECTIVE OPTICAL DEPTH 
DISTRIBUTION FROM Z ~ 5.8-5.4 

In the preceding sections, we attempted to explain observations of 
the variations in Lya forest opacity at z ~ 5.6 in the context of a 
radiation field regulated by a fluctuating mean free path. However, 
another interesting feature of the Lya forest at these redshifts is the 
rapid decline in the average transmission above z ~ 5.5 (Fan et al. 
2006) coincident with increased variations in transmission between 
sightlines (B15). Here we examine whether the fluctuating mean 
free path model described above can account for this evolution. In 
particular, we consider models at a = 5.4 and a = 5.8 to compare 
with the 2 = 5.3-5.5 and z = 5.7-5.9 distributions from B15. 

To construct models at different redshifts that are consistent 
with each other, we follow the same procedure from Section 2 


MNRAS 000, 000-000 (0000) 




8 


F. B. Davies, S. R. Furlanetto 




Figure 6. Maps of the 50 Mpc//i-projected Teff in the Lya forest at 2 : = 5.6 centered on the same slice of the density field shown in Figure 4. The left panel 
shows the Teff map for a uniform ionizing background, where the opacity is correlated with the density field (see the upper left panel in Figure 4). The middle 
panel demonstrates the effect of including a fluctuating ionizing background with a uniform A: the relationship between projected density and Teff reverses 
on large-scales due to the correlation of the density field with the radiation field. The right panel includes the full fluctuating-A ionizing background model, 
greatly increasing the correlation and leading to very high Teff in the center of large-scale voids. 



Figure 7. Top: Cumulative optical depth distribution P{< Teff) computed 
for the uniform background (blue curve) and A = 15 Mpc fluctuating back¬ 
ground models with uniform-A (dashed red) and fluctuating-A (solid red). 
The observed distribution from Becker et al. (2015), including lower limits, 
is shown as the black curve. Bottom: The same curves as the top panel but 
recast as P{> Teff) = 1 — P{< Teff) and shown on a logarithmic scale to 
emphasize the high-Teff behavior. 

to create density and halo fields at 2 : = 5.4 and 2 : = 5.8 us¬ 
ing the same initial conditions. Keeping the ratio between ioniz¬ 
ing and non-ionizing UV luminosity fixed, and for a fixed limiting 
UV magnitude (which is roughly consistent with our halo matching 
procedure), integration of the Bouwens et al. (2015b) luminosity 
functions (interpolated between 2 : ~ 5-6) results in an evolution of 
the ionizing emissivity eion oc (1 + 2 :)“^. As a rough approxima¬ 
tion for the evolution of A at fixed Fhi, we assume the power law 
evolution A oc (1 + 2 :)“^'^ (comoving) measured by Worseck et al. 



Figure 8. Top: Cumulative optical depth distribution P{< Teff) computed 
for the uniform background (dashed blue curve) and fluctuating background 
models (solid color curves). The red, orange, and green curves show the 
fluctuating-A = 15, 22, and 34 Mpc models, respectively. The observed 
distribution from Becker et al. (2015), including lower limits, is shown as 
the black curve. Bottom: The same curves as the top panel but recast as P(> 
Tgff) = 1 — P(< Teff) and shown on a logarithmic scale to emphasize the 
high-Teff behavior. 


(2014b) across 2 < 2 : < 5 where Phi is roughly constant (Becker 
& Bolton 2013). 

Primarily as a result of the rapid evolution in A, in this model 
the average Phi increases by more than a factor of two from 
2 ; = 5.8 to 2 = 5.4. In Figure 10, we show the resulting evolu¬ 
tion of the cumulative Teff distributions in our A = 15 Mpc model. 
The solid curves show the “best-fit” distributions obtained by tun¬ 
ing the n parameter (equation 7) separately at each redshift. The 
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Figure 9. Relationship between average density along simulated sightlines 
and the resulting effective optical depth in the uniform (purple) and fluctuat¬ 
ing (solid: fluctuating-A, dashed: uniform-A) ionizing background models. 
The contour levels enclose 68%, 95%, and 99% of sightlines. When the 
ionizing background is uniform, the lai'ge-scale opacity of the IGM is a 
tracer of the density field, as shown by the tight comelation of the blue con¬ 
tours (see also Figure 6). In contrast, in our fluctuating ionizing background 
model the most opaque regions of the Lya forest correspond to underdense 
regions where the radiation field is suppressed. The addition of mean free 
path fluctuations greatly increases this effect. 

2 = 5.6 and « = 5.8 curves require the same k to fit the oh- 
servations, while a = 5.4 requires a substantially higher value. 
The dotted curve shows the 2 = 5.4 distribution using the same 
K as the other redshifts, demonstrating a clear disagreement with 
the measured values. While the rapid Thi evolution in our model 
is consistent with the observed evolution in the Lya forest trans¬ 
mission from 2 ~ 5.6-5.8, it significantly underestimates Teu at 
2 = 5.4, requiring a substantial adjustment in k from ~ 0.14 to 
~ 0.22 to match observations. The simplest interpretation of this 
K adjustment is that our model overestimates Thi at 2 = 5.4 by 
about a factor of ~ 1.6, essentially requiring Thi to be constant 
from 2 = 5.6-5.4 while increasing sharply from 2 = 5.8-5.6. 


6 DISCUSSION 

Models of the ionizing background that assume either a uniform 
background or a uniform mean free path of ionizing photons are 
unable to reproduce the distribution of Gunn-Peterson troughs 
at 2 >5.4 (B15). Using our semi-numerical ionizing background 
model, we find that introducing fluctuations in the mean free path 
greatly enhances the fluctuations in the strength of the ionizing 
background. The increase in fluctuations manifests on the large 
scales (> A) required by observations. Mean free path fluctua¬ 
tions come about due to the inherently fluctuating ionizing back¬ 
ground in the Universe immediately following reionization - the 
average mean free path is short and early halos are strongly biased, 
so fluctuations in the ionizing background of roughly a factor of 
two above and below the mean are inevitable (Mesinger & Furlan- 



Figure 10. Top: Cumulative optical depth distributions P(< t^s) in th£ 
fluctuating A model (green, red, purple curves) compared to the observa¬ 
tions from Becker et al. (2015) (black curves) at z ~ 5.4, 5.6, 5.8 from left 
to right. The dotted green curve shows the distribution at z = 5.4 without 
renormalizing the optical depths. Bottom: The same curves as above, but 
recast as P(> "Teff) and shown on a logarithmic scale to emphasize the 
high-Teff behavior. 

etto 2009) with underdense regions experiencing a weaker ionizing 
background than overdense regions. 

The resulting mean free path fluctuations then depend on the 
competition between two effects: the regulation of absorbers by the 
ionizing background (MU) and the relative number density of ab¬ 
sorbers due to the large-scale density field. Under reasonable as¬ 
sumptions (A oc r^'^j^A”^), we find that on average the mean free 
path decreases in voids and increases in biased regions, leading to 
enhanced and diminished ionizing background strengths, respec¬ 
tively. Because the radiation field is effectively “smoothed” by the 
mean free path, the effect of mean free path fluctuations is most 
prominent on scales larger than the average mean free path in the 
IGM. 

Applying these enhanced ionizing background fluctuations 
back onto the simulated density field, we find that the typical pic¬ 
ture of the Lya forest where transmission and density are anti¬ 
correlated reverses on large-scales (Figures 6 and 9). The resulting 
distribution of effective optical depths is very similar to observa¬ 
tions if the average mean free path is relatively short (A < 20 Mpc), 
although in detail we find that our model has difficulty matching 
the low-Teff and high-Tefi ends of the distribution simultaneously. 

6.1 Comparison to previous work 

The impact of mean free path fluctuations on the Lya forest has 
also been examined by Pontzen (2014) and Gontcho A Gontcho 
et al. (2014), albeit in the context of baryon acoustic oscillations 
measurements at 2 ~ 2.3. Pontzen (2014) showed that when mean 
free path fluctuations were included the bias of neutral hydrogen 
on scales larger than the mean free path becomes negative (see 
their Figure 2), in qualitative agreement with our results. In con¬ 
trast to our non-linear 3D approach, both of these authors applied 
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linear theory arguments to analytically compute the effect of mean 
free path fluctuations through scale-dependent bias of the radiation 
field. They found that the measured flux power spectrum contains 
valuable information that may allow constraints on the properties 
of sources of ionizing photons and the mean free path (see also 
Pontzen et al. 2014). We leave simulations of the flux power spec¬ 
trum to future work, but note that such arguments may still apply 
at 2 > 5 albeit on much smaller scales. 

As discussed above, our model provides a much better fit to 
observations of the Lya forest at 2 ~ 5.6 than the uniform mean 
free path models in B15. However, a direct comparison is com¬ 
plicated by our different methods for simulating the density field, 
the B15 model having higher resolution and a much more accu¬ 
rate simulation of the gas physics. In contrast to B15, we do not 
find that uniform mean free path models of the fluctuating ionizing 
background produce a wider distribution of Ten than uniform back¬ 
ground models (see Figure 7), but this may be due to our inclusion 
of much fainter sources (Muv ~ 13 vs. Muv ^ — 18 in B15) 
which act to smooth the radiation field because of their weaker level 
of clustering. 

The ionizing emissivity in our fiducial A = 15 Mpc model 
is Cion ~ 1.8 X 10^®((rHi)/10“^^'® s”’^) erg s“^ Mpc“^ Hz“^. 
This compares favorably with the estimates by Becker & Bolton 
(2013) and the required slow evolution to 2 > 6 implied by reion¬ 
ization constraints (Robertson et al. 2015; Bouwens et al. 2015a). 
We cannot claim to constrain the true value of eion with our model 
- in the context of our Lya forest modeling, the low k value we use 
implies a higher Th i by a factor of ~ 7 (as compared to a naive 
extrapolation between the measurements at 2 ~ 5 and 2 ~ 6, e.g. 
Calverley et al. 2011; Wyithe & Bolton 2011), but this ignores the 
substantial error in the computed transmission due to the low res¬ 
olution of our simulated density field (Bolton & Becker 2009) and 
other uncertainties due to the thermal history of the IGM and in¬ 
accuracies inherent to the FGPA method (e.g. Bolton et al. 2005). 
In addition, the average mean free path required to reproduce the 
observed tail to high optical depths may differ depending on the 
F and A dependencies of the mean free path. However, it is reas¬ 
suring that our model is not grossly inconsistent with the current 
understanding of ionizing photon output at 2 ~ 5.5. 

Recently, Chardin et al. (2015) attempted to address the broad 
distribution of effective optical depths at 2 > 5 by computing the 
radiative transfer of ionizing photons through small volume (10- 
40 Mpc on a side) high-resolution cosmological simulations. After 
the end of reionization, the ionizing background in their simulation 
volumes is extremely uniform, resulting in only minor variations in 
Teff between sightlines. As a potential explanation for the broad dis¬ 
tribution observed by B15 they presented a toy “rare source” model 
where the ionizing background was dominated by a population of 
luminous quasars. Such a model is unlikely given the constraints 
on the number density of such sources (e.g. McGreer et al. 2013), 
and we suggest that a much simpler explanation for the uniformity 
of their simulated ionizing background exists. Because their simu¬ 
lation volumes are comparable or smaller than the mean free path 
after reionization (10^0 Mpc), the “smoothing” effect of the mean 
free path on the radiation field naturally leads to an almost com¬ 
pletely uniform background. We believe that high-resolution radia¬ 
tive transfer modeling, similar to that of Chardin et al. (2015), is 
critical to understanding the regulation of neutral gas in the IGM by 
the ionizing background during and after reionization (e.g. Rahmati 
et al. 2013), but modeling the effect of mean free path fluctuations 
(i.e. large-scale coherent variations in the properties of absorbers) 
requires volumes that are at least several A on a side. 



Figure 11. Top: Comparison of cumulative optical depth distributions P(< 
r) for three different realizations of cosmological initial conditions (indi¬ 
cated by colors — red is the fiducial model discussed in previous sections) 
with fluctuating-A = 15 Mpc radiation fields. In this panel, the three real¬ 
izations are nearly identical. Bottom: The same cuives as above, but recast 
as P{> Teff) and shown on a logarithmic scale to emphasize the high-Tefj 
behavior. The different realizations are clearly discrepant at the high-Tefj 
end, indicating that the (400 Mpc)® volume still suffers from cosmic vari¬ 
ance at the factor of ~ 3 level at rgff ~ 7. 


6.2 Cosmic variance 

The volume of the simulations presented in the preceding sections 
is 400 Mpc on a side, corresponding to > 10 mean free paths. 
While this is considerably larger than previous computations of the 
fluctuating ionizing background (e.g. Mesinger & Furlanetto 2009; 
Becker et al. 2015), it is worth examining whether this volume is 
in fact large enough to fully capture the fluctuations seen in obser¬ 
vations. As a simple test, we ran two additional simulations with 
the same model parameters but different realizations of the cosmo¬ 
logical initial conditions. Figure 11 shows the resulting Tefi distri¬ 
butions with fluctuating-A = 15 Mpc radiation fields at 2 = 5.6 
for all three realizations. The distributions are in agreement 
for Teff < 4, but they begin to diverge in the high-Ten tail reaching 
a factor of ~ 3 discrepancy at ~ 7. Note in particular that 
our fiducial model has the weakest tail at high Ten, with the other 
realizations slightly closer to the B15 observations at the highest 
optical depths. 

The importance of cosmic variance in the (rare) high-Tefi tail 
is not surprising. In our model, long wavelength modes of the den¬ 
sity field - reflected in the distribution of sources and the resulting 
radiation field - are the source of Lya forest opacity variations. 
However, there are only ~ (400 Mpc)®/(50 Mpc//i)® ~ 200 
independent regions in our simulated volume at the scale of the 
Teff measurements, limiting the predictive power of each realiza¬ 
tion to features more common than the few-percent level. This ar¬ 
gument ignores the fact that even larger modes also play an impor¬ 
tant role in producing coherent background fluctuations (see, e.g., 
the ~ 200 X 50 Mpc void in Figure 4) which may require yet 
larger volumes to model accurately, beyond the reach of our ex- 
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isting computational resources. Such large features in the radiation 
field may require time-dependent modeling of ionizing source and 
IGM opacity evolution due to their substantial light-crossing time, 
an effect which we ignore in this work. 

6.3 Variation of ionizing source and absorber parameters 

The important assumptions we make for sources of ionizing pho¬ 
tons in our model are as follows: 

(i) Duty cycle of unity 

(ii) Fixed ratio of ionizing to non-ionizing UV luminosity 

(iii) Minimum halo mass of ~ 2 x 10® Mq 

These three assumptions all relate to the effective bias of ioniz¬ 
ing emissivity in our simulations. If the emissivity bias is larger, 
then the resulting fluctuations in the ionizing background and mean 
free path will be stronger. Decreasing the duty cycle from unity 
would result in an increased emissivity bias and thus result in 
stronger fluctuations, although likely only modestly (see: Mesinger 
& Furlanetto 2009). The ratio of ionizing to non-ionizing UV lu¬ 
minosity is typically assumed to be fixed, but the stellar properties 
and environments of UV-faint and UV-luminous galaxies are likely 
to be different. Recent simulations by Wise et al. (2014) found that 
more massive (and thus more luminous) galaxies exhibited smaller 
/esc than lower mass galaxies due to the relative robustness of their 
interstellar medium to supernova explosions. Reducing the ioniz¬ 
ing photon output from high-luminosity galaxies would reduce the 
emissivity bias, leading to weaker fluctuations. However, due to the 
steep faint-end slope of the Bouwens et al. (2015b) UV luminosity 
function, the dominant contributors to the ionizing photon budget 
are already relatively faint galaxies (Muv ~ —18), so the effect 
on our results would likely be modest. For a similar reason, re¬ 
ducing the halo mass cutoff is unlikely to have a substantial effect 
on ionizing background fluctuations because the contribution from 
M\}y > —13 galaxies to the UV luminosity density is small. 

As discussed above, our mean free path model has two key 
parameters: the dependencies on Fhi and A. The mean free path 
should change with the strength of the ionizing background through 
regulation of the size of neutral absorbing clouds. Our assump¬ 
tion of A oc is equivalent to assuming a gas density PDF 

P(A) oc A“® ®. This assumption may be somewhat conservative 
- Ml 1 found a steeper PDF slope at« ~ 5.6 in their hydrodynamic 
cosmological simulations corresponding to A oc By running 

additional simulations with this steeper dependence, we find that it 
has a modest effect that is barely distinguishable from simply ad¬ 
justing the average mean free path. For example, a A = 22 Mpc 
model run with A oc is nearly identical to our fiducial model 
with A = 15 Mpc and A oc 

The mean free path should also depend on the density of the 
local environment - higher density regions will have more dense 
gas to absorb ionizing photons. For simplicity, we assumed that the 
number density of absorbers in a given volume is proportional to 
A, or A oc A“^. The extremely flat evolution of the ionizing back¬ 
ground at 2 ~ 2-5 suggests that neuhal gas is associated with dark 
matter halos (Munoz et al. 2014) so we may be underestimating 
the effect that the density field has on the mean free path. How¬ 
ever, at 2 ~ 6, when the ionizing background is evolving extremely 
rapidly, the overdensities associated with self-shielded gas become 
smaller than the virial overdensity (Munoz et al. 2014), and so the 
bias of absorbers may be small. The sharp change in the inferred 
redshift evolution of the ionizing background at 2 ~ 5.5 (e.g. Fan 


et al. 2006) may reflect a transition from neutral gas residing pre¬ 
dominantly in the IGM to neutral gas residing predominantly in 
collapsed objects. 


7 CONCLUSION 

In this work, we have constructed a 3D semi-numerical model of 
the ionizing background that, for the first time, self-consistently in¬ 
cludes the effect of fluctuations in the mean free path due to the 
ionizing radiation field and density field. We combine the semi- 
numerical halo model of DEXM (MF07) with a prescription for 
mean free path variations motivated by analytic calculations and 
hydrodynamic simulations (Mil) and compute the ionizing back¬ 
ground in a volume 400 Mpc on a side, iterating until the spa¬ 
tially variable mean free path and ionizing background are self- 
consistent. The resulting radiation field shows strongly enhanced 
fluctuations relative to previous models that assume a uniform 
mean free path. Applying the FGPA to sightlines through the quasi- 
linear density field in our simulation volume, we find that the addi¬ 
tion of mean free path fluctuations substantially broadens the dis¬ 
tribution of Lya forest effective optical depths on large scales, par¬ 
ticularly at the high-Teff end which eluded previous work. 

Our results show that the end phases of reionization are rich 
in complexity and astrophysics: there is not a sharp transition from 
the epoch of reionization to the simple uniform background that is 
consistent with the Lya forest at lower redshift. Rather, the inho¬ 
mogeneous source and IGM distributions leave observable artifacts 
even after reionization is over. In principle, we can use these signa¬ 
tures to learn about the absorber population at the tail end of reion¬ 
ization and better understand how they regulate the final phases of 
that process. If our understanding of the IGM improves, this may 
allow us to model the disappearance of the last neutral regions in 
the Universe. The “self-shielding” of void regions by absorbers at 
their edges (e.g. Crociani et al. 2011) may prolong the end stages of 
reionization by acting as large-scale sinks of ionizing photons, sim¬ 
ilar to the semi-analytic modeling by Sobacchi & Mesinger (2014) 
which showed the importance of a F-dependent recombination rate. 

Similar effects may be present in the Hen Lya forest after 
Hen reionization at 2 ~ 3. Indeed, large variations in are seen 
on ~ 40 Mpc scales down to 2 ~ 2.7, and regions with relatively 
low Teft appear to exist as early as 2 ~ 3.5 (Worseck et al. 2014a). 
These observations have been interpreted as an early and extended 
Hen reionization process (Compostella et al. 2014). In contrast to 
the hydrogen ionizing background model in this work, the domi¬ 
nant source of He n ionizing background fluctuations is the rarity of 
bright quasars that likely dominate the He ll-ionizing photon bud¬ 
get (e.g. Furlanetto 2009). These fluctuations may similarly “seed” 
fluctuations in the mean free path, leading to large-scale variations 
in Teff that mimic the end of reionization (Davies, Furlanetto, & 
Dixon, in prep.). 

Many avenues exist to improve the predictive power of our 
model. Potential modifications include a larger simulation volume 
to better characterize the effect of cosmic variance and long wave¬ 
length modes in the density field, high-resolution N-body and/or 
hydrodynamic simulations for a more realistic IGM, and more so¬ 
phisticated models for variations in the mean free path as a func¬ 
tion of environment. Understanding the radiative feedback between 
ionizing sources and neutral absorbers may be critical to explain¬ 
ing the sudden shift in Lya forest evolution at 2 ~ 5.5 (Fan et al. 
2006; Becker et al. 2015; Munoz et al. 2014), and the distribution 
of effective optical depths is a key piece of that puzzle. 
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In the final stages of preparing this paper, we became aware of 
D’Aloisio et al. (2015), which presents an alternative explanation 
for the broad distribution of Ten at z > 5. In their model an ex¬ 
tended reionization process, ending at a ~ 6, leads to strong spatial 
variations in residual heat, which lead to corresponding large-scale 
variations in TeS (tgp oc see equation 7). The resulting 

relationship between overdensity and Ten is reversed compared to 
our model: voids, which are ionized later, remain hot (i.e. trans¬ 
parent), while dense regions, which are ionized first, have had time 
to cool and become more opaque (although this relationship is not 
taken into account in their transmission spectra). This effect should 
in principle dilute the effect of a fluctuating ionizing background, 
and it remains to be seen which effect (if either) dominates. 

Because of the direct connection between large-scale features 
in the galaxy population to Gunn-Peterson troughs in the IGM, 
it may be possible to directly test our model through surveys of 
2 ~ 5.5-5.9 galaxies in z > 6 QSO fields. For example, comparing 
the distribution of galaxies in the upper right panel of Figure 4 to 
the effective optical depth map in Figure 6, the number of > L* 
galaxies in fields corresponding to sightlines with large opaque 
troughs at the same redshift is much smaller than fields with ex¬ 
cess transmission. In future work we will use our model to make 
predictions for the correlation between large-scale Lya forest Ted 
and UV-selected galaxy populations that can be directly tested with 
current ground-based observatories. In principle these predictions 
would allow a definitive test of our model versus the D’Aloisio 
et al. (2015) model. 
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